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SECOND-ORDER  STATISTICS  OF  SPECTRAL  AND  CORRELATION  ESTIMATES 
OBTAINED  BY  MEANS  OF  WEIGHTED  OVERLAPPED  FFT  PROCESSING 


1.  INTRODUCTION 


An  efficient  method  of  estimating  the  spectral  and/or  correlation  characteristics  of  a 
stationary  random  digital  temporal  process,  {g(k)},  is  by  the  use  of  weighted  overlapped  fast 
Fourier  transform  (FFT)  procedures.  In  particular,  a  weighting  sequence  {w(k}}  of  finite 
duration  K  is  overlaid  on  the  long  input  data  sequence  (g(&)}  and  the  product  is  subjected  to  an 
FFT  of  size  N,  yielding  complex  frequency  coefficients.  Parameter  N  governs  the  frequency 
spacing  of  the  coefficients.  The  magnitude-squared  frequency  coefficients  for  this  weighted 
segment  of  data  are  stored  in  computer  memory.  Then,  the  weighting  sequence  is  lagged 
(delayed)  by  L  units  of  time  and  the  FFT  procedure  is  repeated  on  the  weighted  input  data. 

When  a  sufficient  number  of  these  frequency-transformed  segments  are  available,  each  with  a 
different  time  lag  L,  an  average  is  formed  of  the  magnitude-squared  frequency  coefficients, 
thereby  yielding  an  estimate  of  the  power  density  spectrum  of  the  input  data.  If  and  when  an 
estimate  of  the  correlation  function  of  the  input  data  is  also  desired,  an  inverse  FFT  of  the 
magnitude- squared  frequency  coefficients  is  performed,  yielding  correlation  estimates  in  the 
time-delay  domain. 

The  joint  higher-order  probability  density  functions  of  the  spectral  and/or  correlation 
estimates  are  often  of  importance  in  signal  processing.  However,  there  are  features  of  the  FFT 
processing  above  that  make  this  virtually  impossible  analytically.  Even  if  input  data  (g(£)}  are 
zero-mean,  white,  and  Gaussian,  the  overlap  of  adjacent  weighting  sequences,  that  is,  lag  L  less 
than  segment  duration  K,  coupled  with  the  nonlinear  operation  of  magnitude-squared  frequency 
coefficients  and  the  subsequent  averaging,  leads  to  statistical  problems  that  are  analytically 
intractable.  Although  the  joint  moment-generating  function  can  often  be  calculated  for  some  of 
these  types  of  processing  operations,  the  required  matrix  manipulations  have  storage  and  size 
requirements  that  are  frequently  not  practical.  Also,  the  interaction  of  the  fundamental  parameters, 
namely,  weighting  duration  K,  temporal  lag  L,  and  FFT  size  N,  lead  to  complications  in  the 
analysis. 

Accordingly,  the  approach  here  is  limited  to  determination  of  the  first-order  and  second- 
order  statistics  of  the  spectral  and  correlation  estimates.  Several  cases  are  considered,  including 
both  white  and  colored  Gaussian  input  data.  Lag  L  is  arbitrary,  thereby  allowing  for  study  of  the 
adjacent  time  segments,  as  well  as  the  separated  (disjoint)  time  segments  of  data.  The  FFT  size 
N  is  also  arbitrary,  thereby  allowing  for  the  effects  of  wraparound  on  the  estimates  and  their 
stability.  In  particular,  the  dependence  between  estimates  obtained  at  different  segment  locations 
and  at  different  frequency  bins  is  investigated,  and  closed- form  results  for  the  means  and 
covariances  are  obtained. 
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Section  2  of  this  report  contains  the  derivation  of  the  second-order  statistics  of  the 
magnitude-squared  spectral  estimates.  Then,  three  different  cases  of  correlation  estimates  are 
investigated  in  sections  3, 4,  and  5,  each  progressively  more  general  and  more  difficult 
analytically.  Section  6  treats  the  statistics  of  the  complex  frequency  coefficients  themselves. 
Finally,  appendixes  A  and  B  present  detailed  investigations  into  several  special  topics. 


2.  STATISTICS  OF  MAGNITUDE-SQUARED  SPECTRAL  ESTIMATES 


Input  data  sequence  (g(£)}  consists  of  independent,  real,  Gaussian,  random  variables  (RVs) 
with  zero  mean  and  unit  variance.  Real  weighting  sequence  {w(k)}  is  defined  for  all  integer  k, 
but  is  nonzero  only  for  k  =  1 :  K.  A  complex  spectral  estimate  in  frequency  bin  n  is  obtained 
according  to 

z  (ri)  =  ^  g  (k)  w(k )  exp(-i27rnk/N)  for  n  =  0:  N  -l.  (1) 

k 

The  infinite  sum  on  k  is  automatically  terminated  to  K  terms  by  the  weighting;  this  isolates  a 
segment  of  the  data  stream  {g(£)} .  The  integer  N  governs  the  frequency  spacing  of  the  complex 
spectral  estimates  {z(«)}.  Since 

z  (N  -  n)  =  z  *  («)  for  n  =  1 :  N  - 1,  (2) 

it  is  only  necessary  to  consider  the  range  «  =  0 :  V/2  for  {z («)},  where  N  is  presumed  even. 

In  addition,  another  set  of  complex  spectral  estimates  is  obtained  from  a  lagged  data  set 
relative  to  sequence  {g(&)},  according  to 

y («)  =  g (k)  w{k  -  L)  exp[-/2;r  n  (k  -  L)/ N]  =  ^  g (k  +  L)  w(k )  exp(-  iln  nk/N),  (3) 

k  k 

where  L  is  the  amount  of  lag  between  the  two  sets  of  data  and  is  arbitrary.  If  lag  L  is  greater  than 
or  equal  to  weighting  duration  K,  the  two  sets  of  spectral  estimates  in  equations  (1)  and  (3)  are 
statistically  independent  of  each  other,  because  they  involve  different  sets  of  independent, 
Gaussian  RVs.  Interest  here  will  be  concentrated  on  the  case  0  <L<K,  meaning  that  the  two 
segmented  data  sets  encountered  in  equations  (1)  and  (3)  are  statistically  dependent  on  each  other, 
thereby  making  complex  RVs  {z («)}  and  {y («)}  statistically  dependent  on  each  other.  It  is  also 
assumed  that  N>  K;  thus,  the  various  integers  of  the  following  spectral  analysis  will  always 
satisfy 

0  <  L  <  K  <  N;  (4) 

however,  L  can  be  arbitrary  in  general. 

The  magnitude-squared  spectral  estimates  of  interest  are 

q(«)  =  |zO)|2  =  X  8(*)  §0')  w(k )  W(J )  exp [-i2nn(k- j) I N]  for  n  =  0 :  N/2,  (5) 

kj 
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and 


P(«)  =  |y(«)|2  =  X  %(k  +  L)  8 (j  +  L)  ™(k)  w(j)  exp \-i2nn(k-  j)/N]  forn  =  0  :  N/2,  (6) 


where  integer  n  need  not  be  equal  to  n.  The  autocorrelation  of  real  data  sequence  {g(/r) }  is 
E{g(k)g(j)}  =  S(k-j), 

where  E{ }  denotes  an  ensemble  average,  and  the  Kronecker  delta  is  defined  here  as 


m= 


1  for  k  =  0 
0  otherwise 


Then,  from  equation  (5),  the  mean  of  magnitude-squared  spectral  estimate  q(«)  is  immediately 
£{*l(w)}  =  X  w2(k)  for  all  n.  (S 


Of  course,  spectral  estimates  {p(n)}  have  the  same  common  mean  value. 

The  crosscorrelation  of  the  two  spectral  estimates  in  equations  (5)  and  (6),  for  two  different 
(or  equal)  frequency  indices  n  and  n,  is 

CC  =  E{q(n)  p(n)}  =  X  £{g(*0  g(/)  g(*  +  L)  g U  +  L)}  w(k)  w(j )  w(k)  w{j) 

m  (10 

x  exp[-il7r {n(k  -  j)  +  n(k-  for  n  =  0 :  N/2,  n  =  0 :  N/2. 

Since  data  sequence  {g(k)}  is  Gaussian,  the  ensemble  average  in  equation  (10)  can  be  written, 
with  the  help  of  equation  (7),  as  a  sum  of  three  terms,  namely, 

S(k-j)  S(k-l)  +  S(k-k-L)SU-i-L)  +  S(k-l-DS(j-k-L).  (11 

Substitution  of  the  first  product  of  Kronecker  delta  functions  into  equation  (10)  yields  the  first 
component  of  the  crosscorrelation  CC  as 


CC,=Y.w1(k)w1(k)=  2>2(t)  =£{q(n)}£-{p(/i)l. 


Thus,  the  sum  of  the  two  remaining  components  of  CC  will  be  the  covariance  of  RVs  q(«)  and 
P(«)- 
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Substitution  of  the  second  product  of  equation  (11)  into  equation  (10)  yields  the  second 
component  of  CC  as 


CC2  =  w(k)  w(k  - L )  w(j)  w(j - L )  exp[— /2/r(«  +  n)(k-j)/N] 


k,j 


y  w(£)  w(&  -  L)  exp[-z'2;r(n  +  n )  kj TV] 


(13) 


At  this  point,  it  is  useful  to  define  the  complex  ambiguity  function  (CAF)  of  real  weighting 
sequence  (w(£)}  as 

%{m,n )  =  Y  w(k )  w(k-m )  exp(-i27ink/N )  for  all  m,n.  (14) 

k 


Then, 

CC2  =\x(L,n  +  n)\\  (15) 

The  CAF  %{m,  n)  has  period  N  in  variable  n.  Also, 

X{L,n)  =  0  if  L>K,  z(L,N-n)~Z*(L,n).  (16) 

The  third  component  of  CC  is  obtained  by  substituting  the  third  product  of  equation  (11) 
into  equation  (10),  namely, 

CC3  =  Yj  w(k)  M.k  -  L)  w(j )  w(j  -  L )  exp[-/2^(«  -n)(k-  j)/N]  =  \ %{L,  n-n)\2.  (17) 

kj 

Thus,  as  noted  under  equation  (12),  the  covariance  of  interest  is 

cov{q(«),p(n)}  =  \x(L,n-n)\2  +\%(L,n  +  n)\2  for  n  =  0\  N/2,  n  =  0:N/2.  (18) 

Equations  (18)  and  (14)  are  the  main  results.  Lag  L  is  arbitrary;  however,  if  L  is  taken 
larger  than  or  equal  to  K,  the  CAF  in  equation  (14)  will  be  zero,  meaning  that  all  the  spectral 
estimates  {q(«)}  and  {p («)}  are  uncorrelated  with  each  other  for  all  n,n.  As  a  special  case  of 
equation  (18),  by  setting  L  =  0,  there  follows 

cov{q(«),q(«)}  =  |^(0, n-n)\2  +\%(Q, n  +  n)\2  for  n  =  0  :  N/2,  n  =  0 :  A/2,  (19) 

which  is  the  covariance  of  just  the  spectral  estimates  (q(«)}.  As  a  further  special  case, 

var{q(«)}  =  j(0,0)2  +  \%{Q,2rif  for  n  =  0  :  N/2.  (20) 
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Thus,  for  example, 

var{q(0)}  =  var{q(A/2)}  =  2  j(0,0)2.  (21) 

For  0  <L<K,  define 

F(L,  n)  =  y]  w(k)  w(k  -  L)  exp(-i2n  n  k/N )  for;?  =  0  :  A  — 1.  (22) 

k 

This  quantity  can  be  realized  in  MATLAB  as  an  FFT: 

F(L, :)  =  fft(w(L  + 1 :  K).  *  w(l :  K  -  L ),  A).  (23) 

Then,  the  CAF  is  available  according  to 

X(L,  n )  =  F(L,l  +  mod(«,  N ))  for  all  n.  (24) 

Since  integers  n  and  n  are  limited  to  0 :  A/2,  the  mod  function  comes  into  play  when 

n  =  A/2  and  «  =  A/2,  or  n  +  n  =  N.  (25) 


Finally,  the  mean  of  the  spectral  estimates  in  equation  (9)  can  be  expressed  in  terms  of  the 
CAF  as  E{q(n)j  =  j(0,0)  for  all  n.  An  example  MATLAB  program  for  the  covariance  matrix  is 
shown  in  figure  1 . 

NON-WHITE  GAUSSIAN  INPUT  DATA 

If  the  input  data  correlation  in  equation  (7)  is  replaced  by 

E{g(k )  g(y)}  =  C(k  -  j),  C(0)  =  1,  (26) 

which  is  a  real  even  function,  then  the  mean  in  equation  (9)  is  replaced  by 

A{q(n)}  =  ^C(£  -  j)  w(k )  w(j)  exp[-  ilnnik  -  j)/N],  (27) 

kj 

Letting  m  =  k-j,  there  follows 

E{q(n)}  =  ^  C(m)  w(k )  w(k-m)  exp(- i27t nm / A)  =  ^  e\p(-i27rnm/N)  C(m )  (/>w(m),  (28) 

k.m  m 

where 
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clear  all 
close  all  hidden 


K=1024;  %  Length  of  weighting 

L=512;  %  Lag  (shift)  of  adjacent  segment; 

<=  L  <  K 


N=2048;  %  FFT  size;  K  <=  N 

w=hamming(K) ;  %  Weighting  sequence 

%  w=hann (K) ; 

%  w=hann(K+l);  w=w(2:K+l); 

%  w=hann(K+2);  w=w(2;K+l); 


0 


average=sum (w . A2 ) ; 

F=f f t (w (L+l : K) .*w(l;K-L) ,N) ;  %  F  Depends  on  lag  L 

F2=F.*conj (F) ; 
cov2=zeros (N/2+1 ,N/2+l) ; 
for  na=0:N/2 

for  nb=0:N/2 

td=F2 (l+mod(na-nb,N) ) ; 
ts=F2 (l+mod(na+nb,N) ) ; 
cov2 (na+1 , nb+1) =td+ts ; 

end 

end 


Figure  1.  MA  TLAB  Example  for  the  Covariance  Matrix 


<f>w  (m)  =  2  w(k )  w(k  ~  m )  (29> 


is  the  autocorrelation  of  weight  sequence  and  =  0  for  m  >  K. 


On  the  other  hand,  if  N>2K,  and  if  the  input  data  correlation  C(k)  is  expressed  in  terms 
of  the  real  and  even  input  data  spectrum  S(m )  according  to  (see  appendix  A) 


_  [l  for|A:|<A/2] 

C(k)  =  G(k)  exp(/2^ km/N)  S(m),  G(i)=  ,  (30) 

m  [0  for  \k\  >  N/ 2) 


and  if  window  W ( n )  is  defined  as 


W(ri)  =  ^  exp(—  i2jt  nk / A)  w(k),  (31) 

k 

then  substitution  into  equation  (28)  leads  to  mean  value 

£{q(X)}  =  X  s(m)  |ff(«-w)|2.  (32) 
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The  product  G(m)  <f>w(rri)  that  occurs  upon  substitution  into  equation  (28)  is  equal  to  <f>w(m)  for 
all  m  when  N>2K.  On  the  other  hand,  if  only  the  weaker  restriction  N  >K  is  in  effect,  the 
gating  function  G(m )  partially  truncates  <f>w{m),  thereby  obviating  the  simple  relation  (32). 

Relation  (32)  is  recognized  as  a  frequency  domain  convolution  of  the  input  data  spectrum 
with  the  magnitude-squared  window  of  the  weighting  sequence.  Equation  (28)  is  a  correlation 
domain  expression  for  the  mean  of  the  magnitude-squared  spectral  estimate  q(«),  whereas 
equation  (32)  is  a  frequency  domain  expression  for  this  same  quantity.  However,  relation  (32) 
requires  that  N>2K,  whereas  relation  (28)  only  requires  N  >K. 

For  non-white  Gaussian  input  data,  the  crosscorrelation  CC  of  two  different  spectral 
estimates  q(«)  and  p(«)  is  still  given  by  general  relation  (10).  However,  equation  (1 1)  is  now 
replaced  by 

Em)  g O')  g (k  +  L)tU  +  L))  =  C(t-  j )  C(k  -  Jj 

+  C(k-k-L)CU-j-L)  +  C(k-j-L)C(j-k-L).  (3 


Substitution  of  the  first  product  of  equation  (33)  into  equation  (10)  yields  the  first  component  of 
the  crosscorrelation  CC  as 

CC,  =  2  C(k  -  j )  C(k  -  j)  w(k )  w(  j)  w(k )  w(J)  exp[-/2^  {n(k  -  f)  +  n(k  -  f)}/N] 

kjU 


=  2  C(k  -  j )  w(k)  w(j )  exp[-  i2n  n(k  -  j)/N] 

kj 

*YaC(&-j)  W(K)  w(J)  exp[-  i2nn(k  - £)/ N] 

kJ 


(34) 


=  E{q(n)}  E{v{n)}. 


Thus,  the  sum  of  the  two  remaining  components  of  CC  will  be  the  covariance  of  RVs  q (n)  and 
P(«)- 


Substitution  of  the  second  product  of  equation  (33)  into  equation  (10)  yields  the  second 
component  of  CC  as 

CC2  =  Yj  C(k-k-L)C(J-j-L)  w(k )  w(j)  w(k )  w(j)  exp[-/2^  {n{k  -  j)  +  n(k  -  f)}/N ] 

kjkj 

=  ^  C(k  -k  —  L)  w{k )  w(k)  exp[-  i2n{nk  +  nk)/N ]  (35) 

k,k 

•iLCU-j-L)  w(j )  wU)  exp[+  i2n(nj  +  nj)jN]. 

jj 
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In  the  k,k  sum,  let  k  =  k  -  m  to  get 

Z  C(m-L )  w(k)  w(k-m)  exp[-/2^«^/iV  -  i2nn(k  -  m)/N] 

k,m 

=  ^  exp(i2nnm/ N)  C(m-L )  w(£)  w(k  -  m )  exp[-  /2;r(«  +  n)k/N]  (36) 

=  txp{i2nwnfN )  C(w  -  T)  j(m, «  +  «), 


upon  use  of  equation  (14).  A  similar  procedure  applied  to  the  j,j  sum  in  equation  (35)  leads  to 


CC,  = 


Z  QxpiilKnm/ N)  C(m  -  L )  j(m,  n  +  «) 

m 

Z  exp(/2^-«p/A)  C(p)  %(L  +  p,n  +  n) 


(37) 


The  third  component  of  CC  is  obtained  by  substituting  the  third  product  of  equation  (33) 
into  equation  (10).  Then,  a  procedure  similar  to  that  employed  in  equations  (36)  and  (37)  yields 


CC3  = 


Z  exp(-i27rnm/ N)  C(m  -  L)  %(m,  n-n) 

m 

Z  exp(-  i2nnpj N)  C(p )  j(T  +  p,n-n) 


Finally,  as  noted  under  equation  (34)  and  using  the  evenness  of  real  correlation  C(m),  the 
covariance  of  interest  can  be  modified  into  the  two  following  equivalent  forms: 


(38) 


cov{q(«),p(«)}  = 


+ 


+ 


Z  exp(/2;r nm/N)  C{m)  z(L-m,n-n ) 

m 

Z  exp(i2nnm/N)  C(m)  z(L  +  m,n  +  n) 

m 

Z  exp(—i27rnm/N)  C(m)  z(L  +  m,n-nj\ 

m 

Z  exp(/27rw  m/ N)  C(m)  z(L  +  m,n  +  n) 


(39) 


By  setting  L  =  0,  there  follows 
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2 

cov{q(n),q(«)}  =  ^  exp(-i27rnm/N )  C(m)  x(m,n-n) 

m 

2 

+  exp{i27tnml  N)  C(m)  %(m,n  +  n)  . 

m 

Finally, 

2 

var{q(«)}  =  exp(-i27rnm/N)  C(m )  x(m,0 ) 

m 

2 

+  exp(z'2;r  nm/N)  C(m)  %(m, 2  n)  . 

m 


(40) 


(41) 


When  the  data  correlation  C(m)  is  specialized  to  S(m ),  namely,  white  Gaussian  input  noise, 
these  results  reduce  to  equations  (18)  through  (20). 


RELATION  TO  CORRELATION  ESTIMATES  OF  THE  INPUT  DATA 

Correlation  estimates  of  the  input  data  (g(&)}  can  be  obtained  from  the  magnitude-squared 
spectral  estimates  {q(w)}  according  to  an  inverse  FFT: 

1  £=! 

r (m)  =  — >  exp(i2ft nm/ N)  q(n)  for  m  =  0  \  N (42) 
N  ^ 

Reference  to  equations  (2)  and  (5)  reveals  that 

q (N  -n)  =  q (n)  for  n  =  1 :  N  - 1 .  (43) 

Use  of  this  relation  in  equation  (42)  shows  that  r  (m)  is  real  for  all  m  and  that 

r (N  -  m)  =  r (m)  for  m  =  1 :  N  - 1 .  (44) 

Then,  it  is  only  necessary  to  consider  correlation  estimates 

{r(w)|  for  m  =  0 :  N/2.  (45) 

Let  M  =  N/2  + 1  and  define  vectors 
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(46) 


q  =  [q(0)  q(l)-q(7V-l)f, 

e(m)  =  —  [1  exp(i27rm/N)  ■■■  exp{i27r(N -l)m/N}]T  form  =  0:N/2, 

N 

r  =  [r(0)  r(l)  •••  r(N/2)]T , 
and  NxM  matrix 

e  =  [e(0 )  e(l)--e(A/2)].  (47) 

Then,  by  reference  to  equation  (42), 

r  =  eT  q.  (48) 

The  mean  of  correlation-estimate  vector  r  is 


E{ r}  =  E{ q}  (49) 

in  terms  of  the  mean  of  spectral  vector  q,  given  by  equations  (28)  or  (32).  Let  the  covariance 
matrix  of  Nx  1  vector  q  be  denoted  by  Cqq;  this  quantity  is  available  from  equation  (40).  The 
covariance  of  vector  r  is  then 

cov{r}  =  £{[r-£{r}][r-£{r}f}  =  er  Cqq  e  (50) 

by  use  of  equation  (48).  The  matrix  Cqq  is  NxN  and  is  real,  while  matrix  e  is  NxM  and  is 
complex;  see  equations  (46)  and  (47).  Nevertheless,  the  MxM  covariance  matrix  of  r, 
obtained  via  equation  (50),  is  purely  real.  Recall  that  M  =  N/ 2  +  1. 

Another  set  of  correlation  estimates  can  be  obtained  from  the  alternative  spectral  estimates 
in  equation  (6),  namely, 

1  JV-l 

s (m)  =  —  V  exp(/2;r  run/ N)  p(n)  for  m  =  0 :  N/2 .  (51) 

N  p0 

The  MxM  covariance  matrix  between  sets  (45)  and  (51)  is  given  by 

cov{r,s  }  =  eTCqpe,  (52) 

where  covariance  matrix  Cqp  is  available  in  equation  (39),  and  is  NxN  and  real. 

A  shortcut  is  possible  by  observing  that  equations  (42)  through  (44)  can  be  combined  to 
yield 
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(53) 


1  N/  2-1 

r(m)  =  —  [q(0)  +  (-l)m  q(Ar/2)  +  2  V  cos(27tnm/N)  q(«)]  for  m  =  0  :  7V/2, 

N  tr 

which  involves  only  real  quantities.  Define  Mxl  vectors  (M  =  N/ 2  +  1) 


=  —  [1  2cos(27rm/N)---2cos{2nm(N/2-l)/N}  (-l)mf  form  =  0:N/2 , 
A 

q  =  [q(0)  q(l)---q(iV/2)]r,  p=[p(0)  p(l)  - p(iV/2)]r, 
and  MxM  matrix 


(54) 


/  =  [/(0)  m-f(N/ 2)]. 


(55) 


Then, 


r  =  /rq  and  s  =  /r  p. 


(56) 


The  covariance  matrix  of  vectors  r  and  s  then  becomes 


cov{r,s}  =  fT  C~  f. 


(57) 


Covariance  matrix  C~  is  MxM  and  real,  while  matrix  /  is  also  MxM  and  real;  here, 

M  =  N/2  + 1 .  Thus,  equation  (57)  has  two  advantages  relative  to  equation  (52),  namely,  smaller 
matrices  (by  a  factor  of  two)  and  purely  real  quantities.  The  information  required  for 
construction  of  the  covariance  matrix  C~  in  equation  (57)  has  been  presented  in  equation  (39). 

More  explicit  relations  for  these  correlation-estimate  statistics  are  presented  in  sections  3,  4, 
and  5.  Also,  the  statistics  for  the  complex  spectral  estimates  are  presented  in  section  6. 
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3.  STATISTICS  OF  CORRELATION  ESTIMATES 
FOR  WHITE  DATA  AND  N  >  K 


Sequence  (g(&)}  consists  of  independent,  real,  Gaussian,  random  variables  (RVs)  with  zero 
mean  and  unit  variance.  Real  weighting  sequence  {w(k)}  is  defined  for  all  integer  k,  but  is 
nonzero  only  for  k  =  1 :  K.  A  complex  spectral  estimate  in  frequency  bin  n  is  obtained 
according  to 

z (n)  =  y]  g (k)  w(k)  exp (-ilnnk/N)  for«  =  0:V-l.  (58) 

k 


The  infinite  sum  on  k  is  automatically  terminated  to  K  terms  by  the  weighting;  this  isolates  a 
segment  of  the  data  stream  (g(&)} .  The  integer  N  governs  the  frequency  spacing  of  the  complex 
spectral  estimates  {z («)}.  Integer  A  is  presumed  even. 

In  addition,  another  set  of  complex  spectral  estimates  is  obtained  from  a  lagged  data  set 
according  to 

y (n)  =  ^  g (k)  w(k-L)  exp[-i2n:n(k  -  L)/N]  =  ^  g (k  +  L)  w(k )  exp (-ilnnk/N),  (59) 

k  k 

where  L  is  the  amount  of  lag  between  the  two  sets  of  data  and  is  arbitrary.  If  lag  L  is  greater 
than  or  equal  to  weighting  duration  K,  the  two  sets  of  spectral  estimates  in  equations  (58)  and 
(59)  are  statistically  independent  of  each  other,  because  they  involve  different  sets  of 
independent  Gaussian  RVs.  Interest  here  will  be  concentrated  on  the  case  0  <L<K,  meaning 
that  the  two  segmented  data  sets  encountered  in  equations  (58)  and  (59)  are  statistically 
dependent  on  each  other,  thereby  making  complex  RVs  {z («)}  and  {y(«)}  statistically  dependent 
on  each  other.  It  is  also  assumed  that  N  >  K;  thus,  the  various  integers  of  the  following  analysis 
satisfy 

0  <  L  <  K  <  N,  (60) 


but  L  can  be  arbitrary  in  general. 

The  correlation  estimates  of  interest  are  defined  according  to  inverse  FFTs  as 


1 


r (m)  =  —  ^  |z(m)|  exp(/2;r  nm/N )  for  m  =  0 :  N  - 1 

N  „=o 


(61) 


and 
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(62) 


1 


AM 


q (m)  =  —  J'j  |y(«)|  exp (i2nnm/N)  for  m  =  0:N-L 


n=0 


By  using  the  property  from  equation  (58)  that 
z{N -n)  =  z*(n)  for  n  =  1 :  N - 1, 


(63) 


it  can  be  shown  that  all  the  correlation  estimates  {r {m)}  in  equation  (61)  are  real.  It  then  also 
follows  that 


r(N -m)  =  r (m)  for  m  =  1:  N 

Therefore,  it  is  only  necessary  to  consider  the  correlation  estimates 
{r (m)}  for  m  =  0 :  N/2 


(64) 


(65) 


in  the  following.  A  similar  set  of  properties  holds  for  the  alternative  set  of  correlation  estimates 
(q(w)}  in  equation  (62). 


The  autocorrelation  of  the  white  data  sequence  {g(k)}  in  equation  (58)  is 
E{g(k)g(j)}  =  S(k-j), 

where  the  Kronecker  delta  function  is  defined  here  as 


(66) 


m = 


[l  for  k  =  0  | 
[0  otherwise] 


(67) 


Then,  from  equation  (58),  the  mean  of  magnitude-squared  spectral  estimate  z(n)  is  immediately 
E{\z(nf }  =  ^  w2  (k)  for  n  =  0 :  N  - 1 .  (68) 


Therefore,  the  mean  of  correlation  estimate  r (m)  in  equation  (61)  is 


JV-l 


E{r(m)}  =  Yjw2(k)  ~  X  exp(/'2^-«w/iV)  w2{k )  SN(m), 

k  A  „=0  k 


(69) 


where 


SN(m)  = 


\l  for  m  =  0,±N,±  2N, ...] 
0  otherwise  i 


(70) 
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is  the  periodic  Kronecker  delta  function.  Therefore,  use  of  equations  (65)  and  (70)  reveals  that 
equation  (69)  can  be  simplified  to 


X  w2  ( k )  for  m  =  0 
E{r(m)}  =  -  k 

0  for  m  =  1 :  N/2 


Of  course,  correlation  estimates  {q (m)}  in  equation  (62)  have  the  same  mean  values. 

The  next  statistical  quantity  of  interest  is  the  covariance  between  correlation  estimates 
r  (m)  and  q  (m)  in  equations  (61)  and  (62),  where  integers  m  and  m  maybe  equal  or  not.  The 
crosscorrelation  of  these  two  RVs  is 

CC  =  E{r(m)  q (m)}  =  — !y  ^  £{|z(«)|2  |y(w)|2}  exp[/2^(nm  +  nm)/ N],  (72] 

N  n,n= 0 

Denote  the  ensemble  average  in  equation  (72)  as  A.  Then,  from  equations  (58)  and  (59), 

A  =  X  £{g  W  g O')  g (k  +  L )  g U  +  L)}  w(k)  WU)  w(k)  wC/) 

kju  (73; 

x  exp[-z7;r{w(&-y0  + w(&-/)}/.Ar|. 

Denote  the  ensemble  average  in  equation  (73)  as  B.  Then,  using  the  Gaussian  character  of  data 
sequence  (g(k)}  and  equation  (66),  this  ensemble  average  can  be  expressed  as  the  sum  of  three 
terms,  namely, 

B  =  8(k  -  j)  S(k  -j)  +  S(k-k-L)  S(j  -  j_- L)  +  S(k  -  j_- L)  8(j  -  k- L ).  (74 

Substitution  of  the  first  product  of  Kronecker  delta  functions  into  equation  (73)  yields  the 
first  component  of  the  average  A  as 


Al=y>£Jw1(k)w2(k)=  Yjw2(k)  ■ 


Use  of  this  result  in  equation  (72)  yields  the  first  component  of  the  crosscorrelation  CC  as 

(  \2  1  N- 1  f  V 

CCX  =  ^  w2(k)  — -  ]T  exp[/2^-(«w  +  «m)/7V]  =  ^  w2(k)  SN(m)  SN(m) 

V  k  )  N  n,n=0  \  k  ) 


=  E{r(m)}  E{(\(m)}, 
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by  reference  to  equation  (69)  and  the  comment  under  equation  (71).  Thus,  the  sum  of  the  two 
remaining  components  of  CC  in  equation  (72)  will  be  the  covariance  of  RVs  r (m)  and  q(m). 

Substitution  of  the  second  product  of  equation  (74)  into  equation  (73)  yields  the  second 
component  of  average  A  as 


A2=y£j  M*)  w(k  ~  L)  Mj)  w(j  -  L)  exp[-/2^r(w  +  n)(k  -  j)/N]. 


(77) 


Kj 


Use  of  this  result  in  equation  (72)  yields  the  second  component  of  the  crosscorrelation  CC  as 
CC2  =  w(k)  w{k-L)  w(j )  w(j-L) 


Kj 

J  N- 1 


x — -  V  exp [-i27r(n  +  n)(k~j)/N  +  i27r(nm  +  nm)/N] 
N 


=  2  w(<k  “ Z)  w(i)  Mi -L)SN(m-k  +  j)  SN(m-k  +  j ). 

Kj 

The  second  periodic  Kronecker  delta  function  is  nonzero  when 

j  -  k~m  +  pN, 


(78) 


(79) 


where  p  is  an  arbitrary  integer,  positive  or  negative  or  zero.  Then,  there  follows 

CC2  =  w(k)  w(k  -  L)  ^  w(k  ~m  +  pN )  w(k  -m  +  p N  -  L)  SN(m-  m+  p N).  (80) 

k  p 

The  p  N  term  in  the  periodic  SN  function  can  obviously  be  dropped.  Since  both  m  and  m  are 
limited  according  to  equation  (65),  the  subscript  N  can  also  be  dropped. 


Now,  define  the  function 
h(k,  L )  =  w(k)  w(k  -  L ), 

which  is  of  length  K-L.  Then,  equation  (80)  yields 
CC2  =  8{m-m )  h(k,L)  ^  h(k-m  +  pN,L ) 

k  p 

~  5(m-m )  ^  h(k,L )  \h(k  -  m,L)  +  h(k  -  m  +  N ,L)  +  h(k  -  m  -  N ,L)  +  •  ■  ■]. 

k 

Also,  define  the  function 


(81) 


(82) 
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(83) 


H(m,L )  =  ^  h(k,L)  h(k-m,L )  forallX,m. 

k 

Function  H(m,L )  is  even  in  m  and 

H(m,L )  *  0  only  for  \m\  <  K  -  L.  (84) 

Equation  (82)  can  now  be  expressed  as 

CC2  =  S(m  -  m)[H(m,L)  +  H(m  -  N ,L)  +  H(m  +  N , !)  +  •••].  (85) 

But  since  N  >  K  and  0  <  m  <  N  /  2,  the  only  terms  that  need  to  be  kept  are 

CC2  =  S(m-m)  [H(m,L)  +  H(m- N,L)]  form  =  0:N/2,  (86) 

where  some  of  the  terms  may  be  zero. 

Substitution  of  the  third  product  of  equation  (74)  into  equation  (73)  yields  the  third 
component  of  average  A  as 

4=2  w (*)  w(k  ~  L )  wO')  w(j  ~  L)  exp[-/'2^-(«  -  n)(k  -  j)/N].  (87) 

k,j 

The  only  difference  with  equation  (77)  is  in  the  polarity  of  the  n  term.  Use  of  this  result  in 
equation  (72)  yields,  after  simplifications  similar  to  equation  (78),  the  third  component  of  the 
crosscorrelation  CC  as 

CC3  =X  h(k,L)  h(J,L)  SN(m-  k  +  j)  SN(m  +  k  -  j).  (88) 

kj 


Definition  (81)  has  also  been  employed  here.  At  this  point,  the  argument  is  identical  to  that  used 
above  in  equations  (79)  through  (86),  with  the  end  result  that 

CC3  =  S(m)  S(m)  H( 0,  L)  +  8{m  -  N/2)  S(m  -  N/2)  2  H(N/2 ,  L).  (89) 

When  results  (86)  and  (89)  are  added  together,  the  end  result  is  the  desired  covariance,  as 
noted  under  equation  (76);  namely, 

cov{r(m),q(w)}  =  5{m-m )  [H(m,L)  +  H(m- N.L)] 

(90) 

+  S(m)S(m)H(0,L)  +  S(m-N/2)S(m-N/2)2H(N/2,L) 
for  m,m  =  0  :  N/2.  In  more  specific  terms, 
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(91) 


cov{r(0),q(0)}  =  2  H(0,L), 
cov {r(N/ 2), q(N/ 2)}  =  4  H(N/2,L), 
cov{r(m),q(w)}  =  H(m,L)  +  H(m- N,L)  for m  =  1 : iV/2-1, 
cov{r(w),q(w)}  =  0  for  m*m. 

The  function  H(m,L)  is  defined  in  equation  (83)  as  the  autocorrelation  of  sequence  {h(k,L)}. 
This  latter  sequence  is  defined  in  equation  (81)  as  the  product  of  the  original  weighting  {w{k)} 
and  a  delayed  version  by  L  units.  An  example  MATLAB  program  for  equation  (91)  is  given  in 
figure  2. 


Figure  2.  MATLAB  Example  for  Equation  (91) 
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4.  STATISTICS  OF  CORRELATION  ESTIMATES 
FOR  NON-WHITE  DATA  AND  N  >  2K 


The  results  in  the  previous  section  were  obtained  by  means  of  frequency  domain  relations. 
It  is  now  more  convenient,  for  non-white  data,  to  work  solely  in  the  time  domain.  From 
equations  (58),  (61),  and  (70),  there  follows  the  correlation  estimate 


1  AM 

r (m)  =  —  V  exp(/2^  mn / N) 
N  to 


X  g (k)  w(k)  exp(-i2;r  nkj N) 


=  X  gW  w(^)  T7  X  exp[i2^(m  -  k  +  k)  n/N] 

ft  N  ti 

=  X  g(*)  g(£)  vy(^)  SN(m-k  +  k). 

k,k 


(92) 


For  a  given  value  of  k,  the  only  values  of  k  that  can  contribute  are 


k-k-m  and  k  =  k-m  +  N,  (93) 

presuming  that  N>K  and  that  m  is  limited  to  the  range  0 :  N/2,  as  in  equation  (65). 

Equation  (92)  then  simplifies  to 

r(m )  =  X  g(£)  g(£  -  m)  w(k)  w(k  -  m)  +  X  g(£)  g(£  -m  +  N)  w(k )  w(k -m  +  N).  (94) 

k  k 

When  the  product 

p(k)  =  g(k)  w(k)  (95) 

is  defined,  equation  (94)  can  be  written  as 

r(m)  =  X  PW P(* - m)  +  X  P(*) P (k-m  +  N ),  (96) 

k  k 

which  is  the  sum  of  two  autocorrelations  of  p(&),  one  of  which  is  shifted  by  N  units.  Thus, 
wraparound  can  be  present  in  correlation  estimate  r (m)  when  the  condition  N  >  K  is  imposed. 
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In  this  section,  interest  is  directed  at  the  more  stringent  condition 


N>2K  and  m  =  0:N/2, 


(97) 


for  which  the  last  summation  in  equation  (96)  is  zero  because  the  two  p  functions  involved 
cannot  overlap;  that  is,  there  is  no  wraparound  in  the  time-delay  domain.  Thus,  the  case  of 
immediate  interest  is  given  by  the  two  random  variables 

r (m)  =  ^  g(k)  g(k  -  m)  w(k )  w(k  -  m ), 

k 

q(m)  =  X  &(-  +  L)  +  L-m)  w(k)  w(k  -  m ), 

k 


where  both  m  and  m  are  limited  to  the  range  0 :  Nj2. 

For  non-white  input  data  {g(&)},  and  by  use  of  equations  (27)  and  (30),  the  means  of 
correlation  estimates  r (m)  and  q (m)  are  immediately  available  as 


E{r{m)}  =  C(m )  xv(k)  w(k  -m)  =  C(m )  0w(m), 

k 

E{q(m)}  =  ^  C(m)  w(k)  w(k  -  m)  =  C(m)  <pw (m). 


If  the  length  K  of  weighting  (w(&)}  is  large,  correlation  estimates  (98)  will  be  only  slightly 
biased  near  their  origins. 

The  crosscorrelation  of  the  two  correlation  estimates  in  equation  (98)  is 
CC  =  E{r(m)  q  (w)}  =  ]T  E  (g  (k)  g  (k  -  m)  g  (k  +  L)g(k  +  L-  m )} 

k£  (100) 

x  w(k )  w(k  -  m )  w(J{)  w(k  -  m). 

Using  the  Gaussian  character  of  the  data  {g(£)},  the  ensemble  average  is  given  by 

C(m )  C(m)  +  C(k  -  k  —  L)  C(k  -  k-  L-  m  +  m)  +  C(k  -  k-  L  +  m)  C(k  -  k-  L  -  m).  (101) 

Substitution  of  the  first  product  of  equation  (101)  into  equation  (100)  yields  the  first  component 
of  crosscorrelation  CC  as 
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(102) 


CC,  =  C(m)  C(m )  ^  w(k)  w(k  -  m)  w(k )  w(k  -  m) 

k,k 

=  C(/n)  C(w)  <f>w (m)  ^,(m)  =  £{r(/«)}  E{q(m)}, 

by  reference  to  equation  (99).  Therefore,  the  sum  of  the  two  remaining  terms  of  CC  in  equation 

(100)  will  yield  the  covariance  of  correlation  estimates  r (m)  and  q (m). 

The  second  component  of  CC  is  obtained  by  substituting  the  second  product  of  equation 

(101)  into  equation  (100): 

CC2  =  ^  C(k  - k-L )  C(k  -  k-  L  -  m  +  m)  w(k )  w(k  -  m )  w(k)  w(k  -  m) 

k,k 

-  ^  C(p-L)  C(p  -  L-m  +  m)  w(k)  w(k-m )  w(k  -  p)  w(k-  p-m)  (103) 

P  k 

=  C(p-L)  C(p  -L-m  +  m)  6(m,p,m  +  p), 

p 

where  substitution  p  =  k-k  was  made,  and  the  third-order  autocorrelation  of  the  weighting 
{w(k)}  is  defined  as 

6{m ,  n,p)  =  ^  w(k )  w{k  -  m)  w(k  -  n )  w(k  -  p).  (1 04) 

k 

The  approximate  evaluation  of  this  third-order  correlation  is  presented  in  appendix  B. 

The  third  component  of  CC  is  obtained  by  substituting  the  third  product  of  equation  (101) 
into  equation  (100): 

CC3  =  C(k  -  k-L  +  m )  C(k  -  k-  L-  m)  w(k )  w(k  -  m)  w(k )  w(k  -  m ) 

k,k 

=  ^  C(p  - L  +  m )  C(p  - L-m )  ^  w(k )  w(A:  -  m)  w(k  -  p)  w{k  -p-m)  (105) 

P  k 

=  ^  C(p  - L  +  m)  C(p-L-m)  0(m,p,m  +  p)- 

P 

Therefore,  according  to  the  observation  under  equation  (102),  the  final  covariance  of  interest  is 

cov  {r (m),  q (m)}  =  ^  \C(k  -  L)  C(k  -  L  -  m  +  m)  +  C(k  -  L  +  m)  C(k  -  L  -  m)\  6{m,  k,m  +  k). 

(106) 
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As  special  cases, 


cov {r(m),r(m)}  =  ^  [ C(k )  C(k  -  m  +  m)  +  C(k  +  m)  C(k  -m)\  0{m,k,m  +  k)  (107) 

k 


and 


var{r(m)}  =  [C2(k)  +  C(k  +  m )  C(k  -m)\  6(m,k,m  +  k). 

k 


(108) 
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5.  STATISTICS  OF  CORRELATION  ESTIMATES 
FOR  NON- WHITE  DATA  AND  N  >  K 


In  this  section,  the  less  stringent  condition,  N>  K,  is  considered;  the  pertinent  relation  for 
the  correlation  estimate  r (m)  is  then  equation  (94): 


r (m)  =  g (k)  g(k  -  m)  w(k )  w(k  -  m)  +  g(k)  g (k  -m  +  N)  w(k )  w{k  - m  +  N ) 

k  k 

=  rl(m)  +  r2(m), 


(109) 


which  exhibits  time-delay-domain  wraparound  in  the  second  term.  Also,  the  additional 
correlation  estimate  q (m)  now  takes  the  form 

q(w)  =  g  (k  +  L )  g  (k  +  L-m )  w(k )  w{k  -  m ) 

k 

+  X!  g(A  +  ^)  g (k  +  L-m  +  N)  w(k )  w(k-m  +  N)  (110) 

k 

=  q,(w)  +  q2(w), 


which  also  exhibits  wraparound  in  the  second  term. 

The  means  of  these  two  correlation  estimates  are  immediately  available  as 

E{r(m)}  =  C(m )  <j)w(m)  +  C(m  -  N )  0w(m  -  N ), 

£{q(w)}  =  C(w)  </>w (m)  +  C(m-N)  <f>w(m-  N), 

upon  use  of  equations  (27)  and  (30). 

The  crosscorrelation  of  estimates  r (m)  and  q(m)  is  given  by 

CC  =  E{r(m )  q(m)}  =  £{r,(m)  qx{m)}  + E{rx(m)  q 2(m)} 

+  E{r2(m )  q,(w)}  +  E{r2(m)  q 2{m}}  =  P  +  Q  +  R  +  S. 

The  first  term,  P,  is  available  from  equations  (109)  and  (1 10)  as 

P  =  £'{g(A:)  g(&  -  m)  g(&  +  L)  g(k  +  Z,  —  w)}  w(£)  w(k  -  m)  w(&)  w(&  -  m). 


(Ill) 


(112) 


(113) 


Using  the  Gaussian  character  of  data  (g(&)},  the  ensemble  average  is  given  by 


Substitution  of  the  first  product  of  equation  (1 14)  into  equation  (113)  yields  the  first  term  of  P  as 


Px  =  C(m)  Cirri )  </>Jm), 


(115) 


by  a  now  familiar  procedure.  Substitution  of  the  second  and  third  products  of  equation  (114) 
into  equation  (113),  and  use  of  the  change  of  variables  p  =  k-k,  leads  to  the  respective  terms 


P2  =  ^  C(p-L)  C(p-L-m  +  m)  0{m,p,m  +  p), 

p 

P3  =  ^  C(p  -L  +  m)  C(p  - L-m )  6{m, p,m  +  p). 


(116) 


For  the  Q  average  in  equation  (112),  the  ensemble  average  and  the  three  components  of  Q  are 


C(m )  C(m  -N)  +  C(k -k-L)  C(k  -k-L-m  +  m-N) 
+  C(k  -  k-L  +  m-N )  C(k  -k-L-  m ), 


(117) 


Qx  =  C(m)  C(m  -  N )  <f>w(m)  <f>w(m  -  N), 

Q2  =  ^  C(p-L)  C(p  -  L  -  m  +  m-  N)  0(m,p,m  +  p  -  N ), 

p 

Q3  =  ^  C(p  -  L  +  m-  N)  C(p-  L-m )  6(m,p,m  +  p-N). 


(118) 


The  corresponding  results  for  the  R  average  in  equation  (112)  are 

C(m  -  N )  C(m)  +  C(k  -k-L)  C(k -k- L-m  +  m  +  N) 

+  C(k  -  k  —  L  +  m)  C)k  -  k-  L  -  m  +  N ), 


(119) 


Rx  =  C(m  -N)  C{m )  <j)w{m  -  N )  <pw{m ), 

R2  =  ^  C(p-L)  C(p  -  L  -  m  +  m  +  N)  0(m-  N,p,m  +  p), 

p 

i?3  C(p-  L  +  m )  C(p  -  L  -  m  +  N)  0(m-  N,p,m  +  p). 


(120) 


Finally,  the  corresponding  results  for  the  S  average  in  equation  (112)  are 

C(m-N)  C(m-N)+  C(k  -  k-  L)  C(k  -k-L-m  +  m-N) 

+  C(k  -k-L  +  m-N)  C(k  -k  —  L-m  +  N), 


(121) 
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(122) 


St  =  C(m  -  N)  C(m-  N )  <j>w  (m  -  N)  <j)w  (m  -  N), 

S2  =  ^  C(p-L)  C(p  -  L-m  +  m )  0{m-N,p,m  +  p-N), 

p 

53=X  C(p  -  L  +  m-  N)  C(p  -  L-m  +  N )  6{m-N  ,p,m  +  p-N). 


The  sum  of  the  first  terms  of  these  quantities,  namely,  P,  +  (2i  +  ^  +  Sl ,  is  recognized  as  the 
product  of  the  means  in  equation  (111).  Therefore,  the  sum  of  all  the  remaining  terms  in 
equations  (116),  (1 18),  (120),  and  (122)  yields  the  desired  covariance: 

cov{r(«7),q(/w)} 

=  ^  [C(p-L)  C(p  -  L  -  m  +  m)  +  C(p  -  L  +  m)  C(p  -  L  -  m)]  0(m,p,m  +  p) 

p 

+  ]T  [C(p  -  L)  C(p  -L-m  +  m  -N)  +  C(p  -L  +  m  -  N)  C(p-L-m)]  0(m,  p,m  +  p-N ) 

p 

+  ^  [C(p-L)  C(p  -  L  -  m  +  m  +  N)  +  C(p  -  L  +  m)  C(p  -  L  -  m  +  TV)]  0{m  —  N,p,m  +  p) 

p 

+  ^  [C(p-L)  C(p  -  L  -  m  +  m)  +  C(p  -  L  +  m-  N)  C(p-L-m  +  N)]  6{m-  N  ,p,m  +  p-N). 

(123) 


At  this  point,  it  is  convenient  to  define  the  function 
T(L,m,m)  =  [C(p-L)  C(p  —  L  —  m  +  m) 

p 

+  C(p  —  L  +  m)  C(p-  L  —  m )]  6{m,p,m  +  p). 


(124) 


Then,  the  covariance  in  equation  (123)  can  be  expressed  in  a  compact  symmetric  form  as 

cov {r (m),  q {m}}  =  T(L,m,m)  +  T(L,m,m-  N)  +  T(L,m-  N,m)  +  T(L,m-  N,m-  N).  (125) 


An  alternative  (slightly  simpler)  expression  for  function  T(L,m,m)  is 

T(L,m,m)  =  ^  [ C(k )  C(k  -  m  +  m)  +  C(k  +  m)  C(k  -m)]6{m,k  +  L,m  +  k  +  L).  (126) 

k 

As  special  cases  of  the  above,  there  follows 

cov (r (m),  r (m) }  =  T(0,  m,  m)  +  T(0,  m,  m  -  N)  +  T(0,  m  -  N,  m)  +  T(0,  m  -  N,  m  -  N)  (127) 


var{r(m)}  =  T(0,m,m)  +  T(0,m,m-  N)  +  T(0,m-  N,m)  +  T(0,m-  N,m  —  N). 


(128) 


The  function  T(L,m,m )  in  equation  (124)  involves 


9(m, p, m  +  p)  =  ^  w{k)  w{k  -  m)  w(k  -  p)  w(k  -  p-m) 

*  (129) 

=  2_,  h(k,m )  h(k  -  p,m)  =  H_(P->m>m)> 

k 

where  auxiliary  function 

h(k,m)  =  w(k)w(k-m).  (130) 

The  sequence  { h(k,m )}  is  a  product  of  two  weighting  sequences,  one  of  which  is  delayed  by  m 
units.  Thus,  0(m,p,m  +  p)  can  be  interpreted  as  a  crosscorrelation  of  {h(k,m)}  with  {h(k,m)}, 
versus  delay  p.  Then,  the  sum  on p  in  equation  (124)  yields 

T(L,m,m)  =  J^[C(p-L)  C(p-L  -  m  +  m)  +  C(p  -  L  +  m)  C(p-  L-  m)\HJ  p,m,m) 

JL,  r  ,  (131) 

=  ^  [ C(k )  C(k  ~m  +  m)  +  C(k  +  m )  C(k  -  m)\H_(k  +  L,m,m). 

k 

For  white  input  data  {g(&)},  C(k)  =  S(k),  yielding 

T(L,m,m)  =  S(m-m)  H_(L,m,m)  +  S(m  +  m)  H{m  +  L,m,-m).  (132) 

As  a  further  special  case  for  m  =  m  =  0,  there  follows 

H(L,0,0)  =  ^h(k,0)h(k-L,0)  =  ^jw2{k)w2{k-L).  (133) 

k  k 

More  generally, 

H_(L, m.m)  =  h(k, m )  h(k  —  L,m)  =  H(L, m ),  (134) 

k 

while 

T(L,m,m)  =  H(L,m),  except  for  m  =  m-  0.  (135) 
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6.  STATISTICS  OF  COMPLEX  SPECTRAL  ESTIMATES 


Section  2  treated  the  statistics  of  the  magnitude-squared  spectral  estimates.  Here,  the 
quantities  of  interest  are  the  complex  spectral  estimates 

z(n)  =  ^  g(&)  w(k)  exp(-i27rnk/N )  for  n  =  0  :  N -l  (136) 

k 

and  the  lagged  complex  spectral  estimates 

z. (n)  =  ^  g (k  +  L )  w(k)  exp(-/2^ n k/N )  forn  =  0:TV-l  (137) 

k 

for  white,  Gaussian  input  {g(&)} . 

For  later  use,  define  window  functions 

W\  (»)  =  X!  w(k)  w(k  -  L )  exp(-/  2  nnk/N)  for  all  n  (138) 

k 

and 


W2 (n)  -  w2  ( k )  exp(-/ 2  7tnk/N)  for  all  n.  (139) 

k 

Both  windows  have  period  TV  in  n  while  their  real  parts  are  even  in  n  and  their  imaginary  parts 
are  odd  in  n. 

The  mean  of  z (n)  is  zero,  while  the  two  second-order  moments  of  z (n)  are 
E{z2(n)j  =  ^  w2(k)  exp(-i47rnk/N )  =  W2(2ri), 


k 

£{K«)f}  =  I>2(* 0  =  ^(0). 

k 

The  two  second-order  moments  for  separated  frequency  bins  n  and  n  are 
E{z(n)  z («)}  =  w2(k )  exp[-/ 2 7r(n  +  n)k/N]  =  W2 (n  +  n), 

k 

E {z («)  z  * («)}  =  w2  (A:)  exp[-/  2  7t  (n  -  n)  k/N]  =  W2  {n  -  n). 

k 


(140) 


(141) 


These  equations  give  the  complete  covariance  information  about  the  complex  spectral  estimates 

{z (»)}• 


27 


Let  the  complex  spectral  estimate  be  represented  in  terms  of  its  real  and  imaginary  parts 
according  to  z (n)  =  x(n)  +  i  y {n).  Then,  £{x(w)}  =  0,  E{y(n)}  =  0.  Upon  substitution  into 
equations  (141),  and  balancing  of  real  and  imaginary  parts,  there  follows 


E{x(n)  x(«)}  =  ^  W2r  (n  -  n)  +  ^  W2r  (n  +  «), 
E  {y  (n)  y(«)}  =  ^W2r(n-n)-^  W2r  (n  +  «), 
£{x(n)  y («)}  =  -|  W2i  (n-n)  +  ~  W2i  (n  +  n\ 
E {x(n)  y(n)}  =  ^  W2i (n-n)  +  ~  W2i  ( n  +  n). 


(142) 


These  equations  contain  complete  covariance  information  about  the  real  and  imaginary  parts  of 
the  complex  spectral  estimates  (z («)}.  If  n  and  n  are  switched  in  the  fourth  expression,  the 
result  is  equal  to  the  third  expression,  because  W2i  is  an  odd  function.  If  sequence  {W2(n)}  in 
equation  (139)  is  computed  by  means  of  an  Appoint  FFT,  then  replace  n±n  by  mod(«  ±  n,N). 


The  second-order  statistics  of  lagged  estimate  z(n)  in  equation  (137)  are  identical  to  the 
corresponding  terms  presented  above.  Now,  let  n  and  n  be  different  (or  equal)  integers.  Then, 

E{z{n)  z (n)}  =  ^  E{g(k)  g (k  +  L)}  w(k )  w(k )  exp[-/2;r  (nk  +  nk)/N] 

k,k 


=  ^  w(k )  w(k-L)  exp[-i2 7T k/ N  -  i2 n n(k  -  L)/ N] 

k 

=  exp(z  2  ffnL/N)^  w(k )  w(k  -  L )  exp[-/  2  n  (n  +  n)  k/N] 

k 

-  exp(/2 n nLj N)  Wx(n  +  n )  =  A(n,n  +  n). 


(143) 


Also,  in  a  similar  fashion, 

E{z(n)  z.  *(«)}  =  ^  w{k)  w(k-  L)  Qxp[-i2nn k/N  +  i2nn(k  -L)/N\ 

k 

=  exp(-/ 2  ft  nL/ N)  Wx(n  -  n)  =  B(n, n  -  n). 


(144) 


Equations  (143)  and  (144)  contain  complete  covariance  information  about  the  complex  spectral 
estimates  {z(n)}  and  (z(«)}  for  arbitrary  lag  L.  That  is,  temporal  overlap  is  allowed  and 
accounted  for  in  these  relations. 

In  terms  of  their  real  and  imaginary  components,  let 

z(ri)  =  x(n)  +  i  y («),  z{n)  =  u(«)  +  i  v(«).  (145) 
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The  means  of  all  four  real  quantities  are  zero.  Substitution  of  equation  (145)  into  equations 
(143)  and  (144),  and  balancing  of  the  real  and  imaginary  parts,  leads  to  the  second-order 
relations: 


E  (x(»)  u («) }  =  -^  Br  («,  n  -  ri)  +  ^  Ar  (n,  n  +  n), 

E (y(«)  v(«)}  =  ^  Br (n,n-n)~ ^  Ar  («, n  +  n), 
E{x(n)  v(«)}  =  -  Bt(n, n  - n)  +  -j  Ai ( n , n  +  n ), 
E{y(n )  u («)}  =  ^  B,  {n,  n  -  n)  +  ^  At  (n,  n  +  n). 


(146) 


For  this  case,  the  last  two  expressions  are  not  equal.  In  terms  of  the  windows  defined  in 
equations  (138),  (143),  and  (144),  there  follows 

A(n,n)  =  exp(i2trnL/N)  W}(n),  B(n,ri)  =  QX$(-i27rnLlN)Wx{ri).  (147) 

Equation  (146)  contains  complete  covariance  information  about  the  real  and  imaginary  parts  of 
complex  spectral  estimates  {z («)}  and  {z («)},  for  arbitrary  lag  L. 


29  (30  blank) 


7.  SUMMARY 


The  first-order  and  second-order  statistics  for  spectral  estimates  and  correlation  estimates 
obtained  by  means  of  weighted  overlapped  FFT  processing  have  been  obtained  under  a  variety  of 
conditions,  including  colored  input  data  and  wraparound  in  the  time-delay  domain.  These  results 
often  require  the  definitions  of  new  auxiliary  functions  that  involve  the  weighting  sequence 
{w(k)},  its  length  K,  the  amount  of  lag  L  between  data  segments,  and  the  size  N  of  the  FFTs 
involved.  In  particular,  a  third-order  correlation  function  of  the  weighting  sequence  was  required 
to  evaluate  some  of  the  covariances  of  interest.  These  auxiliary  functions  must  then  be  applied 
in  further  summations  to  obtain  the  final  closed- form  results  for  the  desired  covariances. 
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APPENDIX  A 

INVERSE  DISCRETE  FOURIER  TRANSFORM  PROPERTIES 


Function  C(k )  is  nonzero  only  for  |ft|  <  Kc.  Define  discrete  Fourier  transform 

SN  (m)  =  y]  exp(-72;r  mk/N)  C(k )  for  all  m,  (A- 1 ) 

k 

where  the  sum  is  over  all  k.  Function  SN  (m)  has  period  N  in  m. 

Now,  define  the  inverse  discrete  Fourier  transform 

C N  (k)  —  — ~  exp(/2;r  km/N)  SN(m)  for  all  k,  (A-2) 

N  m(N) 

where  the  sum  is  over  any  length  N  interval  (one  period)  in  m.  Function  CN  (k)  has  period  N 
in  k;  in  fact, 

CN  ( k )  =  C(k)  +  C(k  ±  N)  +  C(k  ±  2  N)  +  ■  ■  ■.  (A-3) 

That  is,  CN  ( k )  is  an  aliased  version  of  original  function  C(k).  However,  if 


N>2KC, 


(A-4) 


then  none  of  the  aliased  lobes  in  equation  (A-3)  overlap  the  origin  lobe  C(k).  Then, 


'C,,(k)  for  Ikl  <  Kc 

=  i  A  ,  ,  ’\  =  G(k)C„(k)  for  all*, 

0  for  \k\>Kc 


(A-5) 


where  “gate  function” 

1  for  |A|  <  Ar/2l 


G{k)  = 


0  for  |k|  >  N/ 2j 


(A-6) 


Thus,  if  N  >  2 Kc,  it  is  permissible  to  express 


C(k)  =  G(k)  —  ^  exp(i27r  km/ N)  SN(m)  for  all  k.  (A-7) 

A  m(jv) 

The  gate  function  G(k )  cannot  be  dropped  in  this  relation. 


A-l 


As  a  special  case,  let  m(N )  =  0 :  N  - 1  and  define 


Then,  for  N  >  2 Kc,  equation  (A-7)  becomes 

C(k )  =  G(k )  —  V  exp(/2^"  km/N )  S(m)  for  all  k. 

N  “ 

Again,  the  gate  function  cannot  be  dropped  in  this  relation. 


(A-8) 


(A-9) 


A-2 


APPENDIX  B 

APPROXIMATE  EVALUATION  OF  d(m,n,p) 


The  function  6(m,n,p)  was  defined  in  equation  (104)  according  to 
9{m,  n,  p)  =  7,  w(k )  w(k  -  m)  w(k  -  n )  w(k  -  p). 


(B-l) 


If  weighting  (w(&)}  varies  slowly  with  k,  as  when  sequence  length  K  is  large,  the  sum  on  k 
becomes  approximately 


6(m,  n,  p)  =  Jfifoc  w(x)  w(x  -  m)  w(x  -  n)  w(x  -  p). 

For  example,  if  Hann  weighting 

w{k )  =  y -  y  cos(2 7t k I K)  fork  =  \  .  K, 

then 


(B-2) 


(B-3) 


w(x)  = 


[y-yCos(2 nx/K)  for  0<x</f] 
!  0  otherwise  j 


(B-4) 


Although  the  integral  in  equation  (B-2)  can  be  carried  out  in  closed  form,  it  is  extremely 
lengthy  and  would  be  time-consuming  to  calculate.  An  alternative  approximation  is  to  use 


w(x)  =  exp 


K) 

2 

K2 

'2, 

for  all  x, 


(B-5) 


which  matches  the  Hann  weighting  (B-4)  at  its  peak.  Then,  from  equation  (B-2),  there  follows 


0{m,  n,  p)  =  — ==  exp 


2V/r 

K 


n 


exp 


4  K: 


n 

4 Y1 


(3m2  +3 n2  +3 p2  -  2mn-2np -  2pm) 


{m2  +n2  +  pz  +(m-n)z  +(n-  p)z  +(p-m)2 } 


2  -in 

This  approximation  yields 

0(0, 0,0)  =  — %=  =  0.282  K, 

2  d7T 


(B-6) 


(B-7) 


B-l 


whereas  the  exact  value  is 


0(0, 0,0)  =  =  —  K  =  0.273  K. 

k  128 

Also,  the  integral  of  equation  (B-4)  yields 

[dx  w4  (x)  =  K. 

J  128 

This  suggests  that  the  scale  factor  l/(2^r)  in  equation  (B-6)  be  replaced  by  35/128. 

More  generally,  consider  weighting  function 

w(x)  =  1  -  b  -  b  cos(2 n  xj K)  for  0  <  x  <  K\  w(Kf  2)  =  1 . 

For  b  =  0.5,  w(x)  is  Harm  weighting,  while  for  b  =  0.46,  w(x)  is  Hamming  weighting, 
follows 

w(x) =  exp 

Then,  equation  (B-2)  leads  to 

n2 

-  b  — —  { m 2  +  n2  +  p2  +  (m  -  n)2  +(n-  p )2  +{p-  m)2 
2  K 

In  particular, 

0(0, 0,0)  =  — £=  (=  0.294  K  for  b  =  0.46), 

2  ^2^  b 

whereas  the  exact  value  is 

0(0, 0,0)  =  K  (1-46  +  962  -10Z>3  +4.37564)  (=0.287  /^  for/)  =  0.46). 


9(m,n,  p)  = 


K 


2-J27rb 


exp 


-b 


2  n 

IT2 


f  K) 
l  2 ) 


for  all  x. 


(B-8) 

(B-9) 

(B-10) 

There 

(B-l  1) 

}  -  (B-l 2) 

(B-l  3) 

(B-l  4) 


These  two  functions  of  b  in  equations  (B-l 3)  and  (B-l 4)  are  nearly  equal  for  b  in  the 
neighborhood  of  (0.3,  0.5).  In  fact,  they  are  equal  at  b  =  0.366948. 
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